capture program drop FiltroDatosPoli
program define FiltroDatosPoli 
syntax varlist
	tokenize `varlist'
    local group `1'

drop if abs(totmuj-20)>15

keep if `group' == 1

gen n = 1
collapse (sum) n, by(totmuj)

gen totmuj2 = totmuj^2
gen totmuj3 = totmuj^3
gen totmuj4 = totmuj^4
gen totmuj5 = totmuj^5
gen muj_rd1 = totmuj==19
gen muj_rd2 = totmuj==18
gen muj_rd3 = totmuj==20
gen muj_rd4 = totmuj==21
gen muj_rd5 = totmuj==22

end

********************************************************************************

set seed 1234

use DTA/Datos_ENIA_95_07_Processed.dta, clear

global SubSample all period1 period2 period3 large small sectorA sectorB sectorC sectorD

********************************************************************************
* Residual Bootstrap (18-21):

capture erase DTA/bootstrap_results_21.dta

local N_groups: word count ${SubSample}
matrix coef = J(`N_groups',9,.)

local j = 1
foreach group of varlist $SubSample {

preserve

FiltroDatosPoli `group'
drop muj_rd5

reg n totmuj* muj_rd*
matrix coef[`j',1] = _b[muj_rd1]
matrix coef[`j',2] = _b[muj_rd2]
matrix coef[`j',3] = _b[muj_rd3]
matrix coef[`j',4] = _b[muj_rd4]
predict yhat, xb
predict r, r

save DTA/bootstrap, replace
sleep 1000

qui forvalues x = 1/200 {
use DTA/bootstrap, clear
keep r
bsample
rename r res_`x'
merge 1:1 _n using DTA/bootstrap
drop _merge
save DTA/bootstrap, replace
sleep 1000
}


qui forvalues x = 1/200 {
use DTA/bootstrap, clear
gen y_`x'=yhat+res_`x'
reg y_`x' totmuj totmuj2 totmuj3 totmuj4 totmuj5 muj_rd*
gen b_0 = _b[_cons]
gen b_1 = _b[totmuj]
gen b_2 = _b[totmuj2]
gen b_3 = _b[totmuj3]
gen b_4 = _b[totmuj4]
gen b_5 = _b[totmuj5]
gen pred = b_0+b_1*totmuj+b_2*totmuj2+b_3*totmuj3+b_4*totmuj4+b_5*totmuj5
gen diff_below = (n-pred)*(muj_rd1==1|muj_rd2==1|muj_rd3==1|muj_rd4==1)*(totmuj<20)
gen diff_above = (pred-n)*(muj_rd1==1|muj_rd2==1|muj_rd3==1|muj_rd4==1)*(totmuj>=20)
egen B_`x' = sum(diff_below)
egen M_`x' = sum(diff_above)
drop b_* pred diff_* y_*
save DTA/bootstrap, replace
sleep 1000
}

keep B_* M_*
keep if _n==1
gen i=1
reshape long B_ M_, i(i) j(iteration)
rename B_ B_`group'
rename M_ M_`group'
cap merge 1:1 _n using DTA/bootstrap_results_21 // En la primera ronda no funciona, por eso tiene capture.
cap drop _merge
save DTA/bootstrap_results_21, replace

restore

local j = `j' + 1
di "`group'"
}

preserve
use DTA/bootstrap_results_21, clear
replace i = 21
global vars B_all M_all B_period1 M_period1 B_period2 M_period2 B_period3 M_period3 B_large M_large B_small M_small B_sectorA M_sectorA B_sectorB M_sectorB B_sectorC M_sectorC B_sectorD M_sectorD
foreach i in $vars {
gen `i'_sd   = `i' // Los voy a calcular más adelante al hacer collapse pero necesito dos variables separadas para mean y sd al colapsar.
gen `i'_mean = `i'
}
save DTA/polynomial_21, replace
restore

********************************************************************************
********************************************************************************
* Residual Bootstrap (18-22):

capture erase DTA/bootstrap_results_22.dta

local j = 1

foreach group of varlist $SubSample {

preserve

FiltroDatosPoli `group'

reg n totmuj* muj_rd*
matrix coef[`j',5] = _b[muj_rd1]
matrix coef[`j',6] = _b[muj_rd2]
matrix coef[`j',7] = _b[muj_rd3]
matrix coef[`j',8] = _b[muj_rd4]
matrix coef[`j',9] = _b[muj_rd5]
predict yhat, xb
predict r, r

save DTA/bootstrap, replace

qui forvalues x = 1/200 {
use DTA/bootstrap
keep r
bsample
rename r res_`x'
merge 1:1 _n using DTA/bootstrap
drop _merge
save DTA/bootstrap, replace
sleep 1000
}

qui forvalues x = 1/200 {
use DTA/bootstrap, clear
gen y_`x'=yhat+res_`x'
reg y_`x' totmuj totmuj2 totmuj3 totmuj4 totmuj5 muj_rd*
gen b_0 = _b[_cons]
gen b_1 = _b[totmuj]
gen b_2 = _b[totmuj2]
gen b_3 = _b[totmuj3]
gen b_4 = _b[totmuj4]
gen b_5 = _b[totmuj5]
gen pred = b_0+b_1*totmuj+b_2*totmuj2+b_3*totmuj3+b_4*totmuj4+b_5*totmuj5
gen diff_below = (n-pred)*(muj_rd1==1|muj_rd2==1|muj_rd3==1|muj_rd4==1|muj_rd5==1)*(totmuj<20)
gen diff_above = (pred-n)*(muj_rd1==1|muj_rd2==1|muj_rd3==1|muj_rd4==1|muj_rd5==1)*(totmuj>=20)
egen B_`x' = sum(diff_below)
egen M_`x' = sum(diff_above)
drop b_* pred diff_* y_*
save DTA/bootstrap, replace
sleep 1000
}

keep B_* M_*
keep if _n==1
gen i=1
reshape long B_ M_, i(i) j(iteration)
rename B_ B_`group'
rename M_ M_`group'
cap merge 1:1 _n using DTA/bootstrap_results_22 // En la primera ronda no funciona, por eso tiene capture.
cap drop _merge
save DTA/bootstrap_results_22, replace

restore

local j = `j' + 1

}

preserve
use DTA/bootstrap_results_22, clear
replace i = 22
global vars B_all M_all B_period1 M_period1 B_period2 M_period2 B_period3 M_period3 B_large M_large B_small M_small B_sectorA M_sectorA B_sectorB M_sectorB B_sectorC M_sectorC B_sectorD M_sectorD
foreach i in $vars {
gen `i'_sd=`i'
gen `i'_mean=`i'
}
save DTA/polynomial_22, replace
restore

********************************************************************************

clear
svmat coef, name(Coefs)
gen group = ""
local j = 1
foreach group in $SubSample {
replace group = "`group'" in `j'
local j = `j' + 1
}

egen Bunching21 = rowtotal(Coefs1 Coefs2)
egen Missing21  = rowtotal(Coefs3 Coefs4)
replace Missing21 = - Missing21
egen Bunching22 = rowtotal(Coefs5 Coefs6)
egen Missing22  = rowtotal(Coefs7 Coefs8 Coefs9)
replace Missing22 = - Missing22
drop Coefs*

reshape long Bunching Missing, i(group) j(i)

save DTA/Results_Polynomial, replace

use DTA/polynomial_21, clear
append using DTA/polynomial_22

collapse (mean) *_mean (sd) *_sd, by(i)

foreach group in $SubSample {
rename B_`group'_mean B_mean`group'
rename B_`group'_sd   B_sd`group'
rename M_`group'_mean M_mean`group'
rename M_`group'_sd   M_sd`group'
}

reshape long B_mean B_sd M_mean M_sd, i(i) j(group) string

merge 1:1 group i using DTA/Results_Polynomial
drop _merge

gen stars_B = cond(abs(Bunching/B_sd)>invnormal(1-0.01/2),"***",cond(abs(Bunching/B_sd)>invnormal(1-0.05/2),"**",cond(abs(Bunching/B_sd)>invnormal(1-0.1/2),"*","")))
gen stars_M = cond(abs(Missing/B_sd) >invnormal(1-0.01/2),"***",cond(abs(Missing/B_sd) >invnormal(1-0.05/2),"**",cond(abs(Missing/B_sd) >invnormal(1-0.1/2),"*","")))

save TABLES/Table4_polynomial, replace
//export excel using DTA/polynomial.xls, sheetreplace firstrow(variables)
